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Mapping stationary axisymmetric phase-space distribution 
functions by orbit libraries 
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ABSTRACT 

This is the first of a series of papers dedicated to unveil the mass composition and 
dynamical structure of a sample of flattened early type galaxies in the Coma clus- 
ter. We describe our modifications to the Schwarzschild code of Richstone et al. (in 
preparation). Applying a Voronoi tessellation in the surface of section we are able 
to assign accurate phase-space volumes to individual orbits and to reconstruct the 
full three-integral phase-space distribution function (DF) of any axisymmetric orbit 
library. Two types of tests have been performed to check the accuracy with which DFs 
can be represented by appropriate orbit libraries. First, by mapping DFs of spherical 
7-models and flattened Plummer models on to the library we show that the result- 
ing line-of-sight velocity distributions and internal velocity moments of the library 
match those directly derived from the DF to a precision better than that of present 
day observational errors. Second, by fitting libraries to the projected kinematics of 
the same DFs we show that the distribution function reconstructed from the fitted 
library matches the input DF to a rms of about 15 per cent over a region in phase- 
space covering 90 per cent of the mass of the library. The achieved accuracy allows 
us to implement effective entropy-based regularisation to fit real, noisy and spatially 
incomplete data. 

Key words: stellar dynamics - galaxies: elliptical and lenticular, cD - galaxies: 
kinematics and dynamics — galaxies: structure 



1 INTRODUCTION 



Since the pioneering work of lSchwarzsc hild ( 1979) orbit su- 
perposition techniques have become an important tool in 
the dynamical modelling of spheroidal stellar systems. Sta- 
tionary distribution functions (DFs) of such systems are sub- 
ject to Jeans' theorem and hence depend on the phase-space 
coordinates only via the integrals of motion. In the axisym- 
metric case these integrals are energy E, angular momentum 
along the axis of symmetry L z and for most potentials an ad- 
ditional, non-classical "third integral" 73. Because any set of 
integrals of motion essentially represents an orbit and vice 
versa, the DF can be approximated by the sum of single- 
orbit DFs, with the only adjustable parameter being the 
total amount of light carried by the orbit. The main task 
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that remains to adequately describe hot stellar systems is to 
find an appropriate set of orbits. 

Orbit superposition techniques have be en used to model 
such systems in various symmetries (e.g. |Rjx_et_alJ M971; 
Romaiiowsl^^T^ochanek 1997|: [van der Ma rel et alJ1 1998: 
ICretton et alJ 1 19991: JCap j^llari_et_alj 120021 : IVerolme et alJ 



bool iGebhardt et al J I2003J : Ivan de Ven etal.1 120031) . with 
the goal of determining different dynamical parameters like 
central black hole masses, internal velocity anisotropy or 
global mass-to-light ratios. An orbit library tracing the 
phase-space structure of a trial potential is fitted to ob- 
served photometry and kinematics, to decide, whether or 
not it gives a valid model of the corresponding galaxy. 

In the spherical case there exists a well-known mass- 
anisotropy degeneracy permitting in general convincing fits 
to the projected velocity dis persion a even if the trial poten- 
tial differs from the true one JBinnev fc Mamorj|l982n . With 
complete knowledge of the full line-of-sight velocity distri- 
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butions (LOSVDs) however, it is po ssible to reconstruct 
the D F, given the potential is known IIDeionghe fc Merrittl 
1992). Furthermore, even for th e realistic case where th e po- 
tential is not k nown in advance iMerritt fc S aha (1993) an d 
Gerhard ( 1993) have shown how the information contained 
in the LOSVDs can constrain both, the potential and the 
DF. 

Likewise, in the axisymmetric case. lDehnen fc Gerh ard 

lll993l) have calculated realistic smooth DFs and have shown 
that a similarity close relationship exists between the poten- 
tial and internal kinematics on the one side and the pro- 
jected kinematics on the other side. However, fits of axisym- 
metric li braries still pose some additi o nal o pen questions. 
Recently, IValluri. Merritt fc Emselleml J2004ft discussed the 
indeterminacy of the reconstruction of the potential in gen- 
eral axisymmetric systems from two- or three-dimensional 
data sets, by studying the shap e of the y 2 -contours describ - 
ing the quality of the orbital fit. lCretton fc Emsell em ( 2004) 
argued that even in case of a mathematically non-degenerate 
f(E, L z )-system an artificial degeneracy occurs, caused by 
the discreteness of the orbit library and emphasized the role 
of appropr iate smoothing , altho ugh not providing a definite 
solution. iRichstone et al.l B2004I) critically analysed their ar- 
guments and emphasized that both high quality compre- 
hensive data sets and orbit libraries are needed to achieve a 
reliable modeling of axisymmetric systems. 

In view of this discussion concerning orbit-based dy- 
namical models it seems worthwhile to step back and in- 
vestigate how well orbit libraries represent the phase-space 
structure of a given dynamical system. This includes the 
examination of the choice of orbits, which in the generic 
axisymmetric case is difficult, since part of the phase-space 
structure is unknown due to our ignorance related to J3. A 
key tool for such an analysis are the orbital phase-volumes 
which accomplish the transformation from the relative con- 
tributions of individual orbits to the library, the so called 
orbital occupation numbers or orbital weights, into phase- 
space densities (and vice versa). The availability of such 
phase- volumes actually offers several applications: 

(i) Accurate phase-volumes allow to calculate internal 
and projected properties like density and velocity profiles, 
line-of-sight velocity distributions (LOSVDs) etc. of general 
axisymmetric DFs f(E,L z ,Iz) via orbit libraries. Besides 
the possibility to systematically study the structure of gen- 
eral axisymmetric systems, these profiles provide a direct 
check of the choice of orbits via their comparison with those 
calculated from directly integrating the DFs. 

(ii) From any fitted library one can reconstruct the cor- 
responding DF via the phase-volumes and thus reconstruct 
the DF from any observed early-type galaxy in the axisym- 
metric approximation. 

(iii) If the library is fitted to some reference data con- 
structed from a DF, then |(ii)| allows to investigate how 
closely the fit matches the input DF and to implement an 
effective regularisation scheme, permitting to fit real (noisy) 
data sets. 



components integrated ab out the unknown integrals they 
have b een exploited bv e.g. JRix et al| (119971) , ICretton et alJ 
J1999D and lVerolme fc de Zeeuwl (12002ft . 

The aim of the present paper is to introduce a gen- 
eral implementation for the calculation of individual orbital 
phase- volumes in any axisymmetric potential and, by follow- 
ing applications | (i) | and | (iii) | to prove that our libraries accu- 
rately map given dynamical systems. This directly justifies 
our setup of the library and sets the basis for our project to 
recover the dynamical structure and mass composition of a 
sample of flattened early- type galaxies in the Coma cluster. 
In a subsequent paper we will focus on the question of how 
much smoothing has to be applied in order to get an optimal 
estimate of the dynamical system underlying a given set of 
noisy and spatially inco mplete observational data. The full 
analysis of the data set iMehlert et alj|2000t IWegner et alJ 
2002) will be addressed in a future publication. 

The paper is organized as follows. In Sec. |5]we define 
all quantities related to the library used in the subsequent 
Sections and describe our orbit sampling. Sec. [3] outlines 
the relation between orbital weights and orbital phase-space 
densities. Sec. ^contains a description of our implementation 
to calculate individual orbital phase-space volumes. In a first 
application we calculate internal and projected properties 
of given DFs using orbit libraries in Sec.|^| Sec. |S| discusses 
how the library is fitted to given data sets and in Sec.|7|we 
reconstruct reference DFs from their projected kinematics. 
Finally, Sec.[S]summarizes the results. 



2 THE ORBIT LIBRARY 

Our method to set up the orbit libraries used for the dy- 
namical modelling is based on the procedure presented in 
Richstone et al. (in preparation). There, the reader finds a 
description of the basic properties of the program. In this 
section we define quantities that are used later on in this 
paper. 

In the following we assume that the luminosity density 
v is known. In the analysis of real data it has to be obtained 
by deprojection of the measured photometry. With the stel- 
lar mass-to- light ratio T — (y-) the mass density pi of the 
luminous material follows from v as pi —T v. 

The total mass density p possibly includes a dark com- 
ponent poM and reads 



p — p (r, r c , v c ) = r v + p D M 



(i) 



IVandervoord 1(1984) touched the problem by writing 
down the transformation from cells of integrals to the corre- 
sponding phase-space volumes. However, the resulting rela- 
tions are only suitable for explicitly known integrals, e.g. for 
single orbits only in the rare case all integrals are known. For 



Once the mass-profile is fixed, the potential $ follows 
by integrating Poisson's equation. With $ known, a large set 
of orbits is calculated, sampling homogeneously the phase- 
space connected with $. 

2.1 Spatial and velocity binning 

As described in Richstone et al. (in preparation) we divide 
the meridional plane into bins, equally spaced in sin?? 1 , lin- 
ear in r near the inner boundary r m i n of the library and 

1 Throughout the paper, we use spherical coordinates (r, $,<p), 
with »9 = 0° corresponding to the equatorial plane. If not stated 
otherwise, we use super- or subscripts h,i,j,k as indices, l,m,n 
as exponents. 
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logarithmic at the outer boundary r max , (if not stated oth- 
erwise we use N r = 20 radial bins, N$ — 5 angular bins.) 
For the projection of the library we use the same binning 
as for the meridional plane. Every spatial bin in the plane 
of the sky is subdivided into JV ve i bins linearily spaced in 
projected velocities between — ii max and t> ma x, leading to a 
binsize for the LOSVDs of 



Aulosvd = 2 



iVv. 



(2) 



Even if the potential is spherical, then our spatial binning 
tags an axis of symmetry. Later, when referring to a "minor- 
axis" we always mean the axis i? = 90° of the library. 

2.2 Orbital properties 

Luminosity. For every orbit i in the library we store its 
normalized contribution to the luminosity dL^ in spatial bin 
1 ^ j *s N r x N$, which equals the fraction of time the 
orbit spends in bin j. Let At k denote the fc-th timestep in 
the integration of orbit i, so that 






is the total time elapsed until timestep k and 
J> ' = {k : (r(*?),tf(t? )) G bin j} 



(3) 



(4) 



contains all timesteps during which orbit i is located in spa- 
tial bin j. Accordingly, we can write 



dU = 






(•>) 



with Ti = ^2 Ati being the total integration time of orbit i. 
Given the orbit's weight Wi to the whole library - the 
integrated luminosity along the orbit - the total luminosity 
of the library in spatial bin j reads 



dv = J2 w > dL i- 



(6) 



Internal velocity moments. To obtain the internal ve- 
locity moments (v l r v™vfy of the orbit library, we store for 
each orbit i and each time step Atf the product of veloc- 
ities v l r v™v™ and fractional time At*/T;. All contributions 
in spatial bin j are summed up to yield 



/ l m n\ J 

\V r Vi} v v ) 



EA-f 
v r v# v v 



T, 



(7) 



Thus for the whole library the velocity moments in spa- 
tial bin j follow as 



/ l m n\i 1 V~^ / I m n\i 



(8) 



Projected kinematics. For the projected kinematics of 
the library we record the normalized contribution to the 
kinematics LOSVCH at projected position j and projected 
velocity 1 $J k $C AT ve i for every orbit. Again for the whole 
library the LOSVD reads 



LOSVD 



ik 






w t LOSVD 



,J k 



(9) 



By fitting a Gauss-Hermite series to th e LOSVD - ,fe we 
obtain the Gauss-Hermite parameters llGerhardl 1 19931 : 
Ivan der Marel fc Franxlll993h 



GHP Jfc = {^ k ,v 3 \a ]k ,Hi k ,Hf} 
of the LOSVD. 



(10) 



2.3 Choice of orbits 

To obtain a reliable representation of phase-space it is im- 
portant that any allowed combination of the integrals of 
motion (E, L z ,Is) is represented to some degree of approxi- 
mation by an orbit in the library. The absence of some orbit 
family in the library might misleadingly emphasize certain 
dynamical configurations in the final fit. 

Sampling E and L z . Richstone et al. (in preparation) 
adjust the orbit sampling in (E, L z )-space according to their 
spatial binning. From the requirement that every pair of grid 
bins Ti ^ Tj in the equatorial plane should be connected by 
at least one equatorial orbit with r pcr i = r; and r apo = 
Tj they uniquely derive a grid of orbital energies E and z- 
angular momenta L z . We experimented with doubling the 
number of peri- and/or apocenters per radial bin but found 
that the above described method yields a sufficiently dense 
sampling of the (E, L z )-plane. 

Sampling 1$. It is common practice among the various 
existing Schwarzschild codes to sample ^3 by dropping or- 
bits at given energy E and angular momentum L z from the 
zero- velocity-curve (ZVC, defined by E — L 2 z /{2r 2 cos 2 1?) + 
$(r, 1?)). Richstone et al. (in preparation) use the intersec- 
tions of the angular rays of the meridional grid with the ZVC 
as starting points. This sampling ensures that each sequence 
of orbits with common E and L z contains at least one orbit 
that is roughly confined to the region between the equatorial 
plane and each angular ray of the meridional grid. 

If we consider only potentials symmetrical about the 
equatorial plane with d$/dz > 0, then every orbit even- 
tually crosses the equatorial plane and leaves a footstep in 
the surface of section (SOS) given by the radii r and ra- 
dial velocities v r of the upward equatorial crossings. Orbits 
respecting a third integral show up in the SOS as nested 
invariant curves, sometime s with embedded resonances (e.g. 
iBinnev fc Trem aine 1987). Fig. shows an example of a 
SOS. The dots mark representative points of invariant curves 
obtained by numerically following orbits with common E 
and L z in a flattened Hernquist potential with total mass 
M — 10 11 Mq, scaling radius r B — lOkpc and a flattening of 
q = 0.5 (see Sect. l5.l1 below for further details). 

The SOS encompasses all available orbital shapes and 
a representative sampling of orbits should end up with the 
SOSs homogeneously filled with orbital imprints. Unfortu- 
nately, we are not aware of any simple relationship between 
the drop-point of an orbit on the ZVC and its corresponding 
appearance in the SOS, as long as ^3 is not known explicitly. 
In order to guarantee a representative collection of orbits in 
any case, we sample the orbits as follows. 

In a first step we drop orbits from the (outer) intersec- 
tions of the angular rays of our spatial grid with the ZVC 
as described in Richstone et al. (in preparation). Then, for 
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any pair (E,L Z ) included in the library we choose Nl radii 
n, 1 ^ I ^ Nl equally spaced in log(r) on the equatorial 
plane between r-p Cr i and r apo of the equatorial radial orbit 
with energy E and angular momentum L z . We start with 
the smallest of these radii n and launch an orbit i from the 
equatorial plane with the maximal radial velocity 



2(E-*(n)) 



a 



V mxK {E,L z ,ri) 



(11) 



For the subsequent orbits i' we stepwise decrease v ri i by 
Av r> i' (see equation l|13fl below) until we reach v r ,i' = and 
pass over to the next radius rj+i. 

With (E,L Z ) and (ri,v r ,i) fixed the orbital v#,i is de- 
termined by 



v#,i(E,L z ,ri,v r ,i) = \ 2(E - $(n)) 



r 2 

=J- (12) 



When « r ,i = 0, then v# ti (E,L z ,r t ,v r ,i) — Vmax(E,Lz,ri). 
For each velocity pair we launch an orbit from the equatorial 
plane at the actual r; with the actual velocities v r ,% and iv,i- 
This procedure is repeated for each of the Nl radii. If, at 
a specific launch position, we find an imprint in the SOS of 
a previously integrated orbit which differs from the current 
launch position by less than 10 per cent in radius and radial 
velocity, then we regard this part of the SOS as already 
sampled and discard the orbit. 

The velocity stepsize Av T ,i is set as 



Av r ,i — min{Au L osvD , £ v m ,i-i] ■ 



(13) 



where Ai>losvd is the width of the LOSVD-bins (cf. equa- 
tion (J2J) and 



v m ,i = max {vi : (rf,vf) £ SOS} . 



(14) 



Here SOS denotes the set of the N sos orbital imprints in the 
SOS and i is the index of the actual orbit. We usually take 
£ = 1/3.5. Trying different values for Nl we found Nl — 30 
sufficient to yield a dense filling of the SOS with approxi- 
mately one invariant curve crossing the r-axis of the SOS in 
each of the equatorial meridional grid bins. The correspond- 
ing orbits have large tf-motion and need to be included in 
the library to avoid a radially biased collection of orbits. 

The velocity stepsize is largest for the radial orbits and 
gradually decreases when the SOS is filled with orbits (note 
that v m ,i-i is the maximum of the radial velocities in the 
SOS of the "precursor" orbit i— 1). For the shell orbits, the 
stepsize becomes smallest. The adjustment of the stepsize 
in each step ensures that we sample the more radial orbits 
with a resolution that corresponds at least to the width of 
the LOSVD bins and that the sampling is refined for the 
shell orbits. 

After the above described sampling is done, we measure 



the maximum f s of all r m i n: i/r n 



with 



rmin.i = mm 



{r* : (r!,v!) £ SOS} 



(15) 



and r max ,i defined analogously. To ensure that the sequence 
contains all orbits up to (approximately) the thin shell orbit, 
we complete the library if necessary by launching orbits from 
the equatorial plane with v r = at 
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r [kpc] 

Figure 1. Example of a surface of section for a flattened Hern- 
quist model (details in the text). All orbits have been integrated 
for iVsos = 80 intersections with the SOS. 



where f 3 = r min y/r max y, until f s > 0.9. 

Fig. Q illustrates the dense coverage of the SOS with 
invariant curves after all orbits are integrated for a flattened 
Hernquist potential. 



2.4 Use of the library 

If the relative contribution of each orbit to the whole library, 
the orbital weight w%, is specified, then according to equa- 
tions 0, |JS| and Q the library provides a specific model 
including the LOSVDs, internal density distribution, inter- 
nal velocity moments and so on of this particular orbit su- 
perposition. 

If the library is constructed to test whether or not a 
given trial mass distribution leads to a consistent model 
of an observed galaxy, then the model, in particular the 
LOSVDs, have to be compared to the observations. If the 
comparison turns out not to yield a satisfactory fit, then the 
weights can either be recalculated (see Sec. |S] for details) 
or the actual mass distribution has to be rejected. If on the 
other hand, the fit shows that the actual set of weights seems 
to be a valid model of the galaxy, then one can reconstruct 
the internal velocity structure and DF from the iWj's. 

Conversely, if one has a DF at hand and wants to calcu- 
late e.g. its projected kinematics without going through the 
appropriate integrals, one can assign the orbital weights ac- 
cording to the DF (see Sec.|^J without any fitting procedure 
and analyse the output of the library. This can be usefull 
to study systematically the projected properties of general 
axisymmetric distribution functions depending on all three 
integrals (E,L z ,Is). 

In the following we will make use of both applications 
with the goal to investigate the accuracy of our orbit li- 
braries. 



3 ORBITAL WEIGHTS AND PHASE-SPACE 
DENSITIES 

In order to reconstruct the DF from the library or to cal- 
culate spatial profiles of internal or projected properties of 
some given DF, one needs to convert orbital weights into 
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phase-space densities and vice versa. This section summa- 
rizes the connection between orbital weights and orbital 
phase-space densities under the regime of Jeans' theorem. 



3.1 Phase-space densities of orbits 

Consider a system in which the orbits respect n integrals of 
motion Ii, . . . ,I n . Because the phase-space density of sta- 
tionary systems is constant along individual orbits (Jeans 
theorem) the phase-space density along orbit i is given as 
the orbital weight w% divided by the phase-space volume 
Vi. More formally, let I denote the n-dimensional set of or- 
bital integrals (Ji,. . . ,J„), let V denote the 6-dimensional 
phase-space, V(V) its power set and let £ : T — > V(V) map 
a n-tuple of orbital integrals (ii, . . . ,/„) £ I on to the hy- 
persurface £(ii , . . . , I n ) Q V in phase-space covered by the 
corresponding orbit, 

C(/i, . . . , /») = {p € V : hip) =I U ..., Up) = /„}. (17) 

With Ui C I being the small cell in integral space repre- 
sented by the orbit i, 

Ui = {(Ii,...,I„)&I:h e [Il,i-AIi,i,h,i + AIi,i], 

...,!„ e [In.i - AI n ,i, /„,, + AJ„,i]} (18) 

we define the characteristic function 



4 ORBITAL PHASE- VOLUMES 



Xi 



1 : {r,&,ip, v r ,v#,v v ) e d 
: (r,&,<p,v r ,v.8,v< f ,) £ d 



of the image set 

O t = (J W 

we5(Wi) 



(19) 



(20) 



of Ui in phase-space. The volume of the phase-space region 
represented by orbit i then follows as 



Vi= I X id 3 rd 3 v 



(21) 



and accordingly the phase-space density along the orbit 
reads 



fi 



Vi 



(22) 



3.2 Orbital weights from DFs 

If we reverse the application of equation (1221 , and assign the 
orbital weights according to some given DF /, 



fiVi, 



(23) 



with fi = f(Ii,i, . • . , In,i) now being the DF / evaluated at 
the orbit's position in integral space, then the DF fut of the 
entire library, which consists of the combined contributions 
of all orbits 



/lib = J2 /• 



x, 



(24) 



will be the mapped version of / on to the library. Equa- 
tion 12311 together with equations JfjJ, (jjj and © can be 
used to calculate the LOSVDs, internal velocity profiles and 
density distribution of any axisymmetric DF with known 
potential. 



Two degrees of freedom. iBinnev. Gerhard fc Hut (1985) 
have shown, that for autonomous Hamiltonian systems with 
two degrees of freedom the phase-volume of any orbit can 
be derived from the SOS by integrating the times between 
successive orbital visits of the SOS 



V ^ AE 



T(r,v r ) drdv r , 



(25) 



where T(r, v r ) is the time the orbit needs from (r, v r ) to the 
next intersection with the SOS, and AE defines a small but 
finite cell around the orbit's actual energy E characterizing 
the hypersurface in phase-space represented by the orbit. 

Axisymmetric case. Richstone et al. (in preparation) 
carry over this result to axisymmetric systems and approx- 
imate the phase-volumes as 



V w AL Z AE 



T(r,v r ) drdv r . 



(26) 



SOS 



Here AL Z and AE stand for the range of energies and angu- 
lar momenta represented by the orbit under consideration. 
Equation 126H is valid independent from the orbit beeing 
regular or chaotic. 

Calculating the SOS-integral. In what follows we de- 
scribe our novel implementation of equation I2fcill that im- 
proves on the method of Richstone et al. (in preparation) to 
deliver higher precision phase-space volumes. 

For all orbits in a sequence with common energy E and 
angular momentum L z we obtain a representative sample S 
of the SOS by storing iV sos imprints of each orbit in the SOS 
given by the radial positions and velocities 2 at the times t { 
of the orbital equatorial crossings, 

5 = {(r!,v?) : rt = r(i t fc(s) ),^ = |«f(*? W )|, 

Ei = E,L Sti = L x ,1^8<:N aoa }. (27) 

Typically, we integrate each orbit up to JV sos = 80 inter- 
sections with the SOS and choose N^ os = 60 points for the 
calculation of the phase- volumes randomly out of the whole 
set of intersections. We also store the time intervals 



t(rl,vl) 



,k(s+l) _ ,fc(s) 



(28) 



between two successive intersections. 

Inspection of Fig. shows that only a tessellation ap- 
proach can numerically integrate equation 1261 in the gen- 
eral case including regular, resonant and chaotic orbits. To 
this purpose we decided to perform a Vor onoi tessellation of 
5 using the software of IShewc huk (1996). This tessellation 
uniquely allocates a polygon to each element of S. The edges 
of the polygon are located on the perpendicular bisections 
of pairs containing the element under consideration and one 
of its neighbours and are equidistant to the actual pair and 
a third element. For almost all elements the polygons are 
closed and encompass an area containing the actual element 



2 To reduce the computational effort we take the absolute value 
of the radial velocities, thereby exploiting the symmetry of the 
SOS with respect to the r-axis in our application. 
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and all points that are closer to it than to any other ele- 
ment. The areas enclosed by the polygons completely cover 
the space between the elements and therefore characterize 
the fractional area inside the SOS occupied by each orbit. 

Fig- HI shows the same SOS as Fig. The open circles 
display r and v r at the orbital equatorial crossings. The 
thin lines around these circles mark the Voronoi cells, e.g. 
the polygons allocated to the elements of S and boundary 
points (see below). 

With AA1 beeing the surface area enclosed by the poly- 
gon around (rf , vf) G S the integral expression in the phase- 
volume of orbit i (cf. equation 1261 can be approximated 3 



/ T(r, v r ) dr dv r « V" t{r 
Jsos 



AAf 



(29) 



At the boundary of the distribution of sampled points, 
there may not be enough neighbours around a given ele- 
ment of S to close its polygon. In order to ensure that every 
Voronoi polygon is closed and confined to an area enclosed 
by the ZVC of the SOS (given by v r of equation Hil l we 
construct an envelope around the distribution of sampled 
orbital intersections. In Fig.|5]the points of the envelope are 
marked by the solid dots. They are constructed as follows. 

In a first step, we determine the maximum vo of radial 
velocities in 5 



Wo 



max {v r : (r, v r ) G 5} . 



(30) 



To ensure that no Voronoi cell exceeds below the axis v r — 0, 
all imprints in the SOS with t>f ^ eiio are mirrored about 
the axis v r — (typically e — 0.1). For the rest of the SOS 
we construct an envelope in an iterative loop starting from 



(f ,Vo) = (r,v r ) G S, v r 



Vo- 



(31) 



In each iteration n+1 we search for (r n +i, V n +i) G 5 obeying 

v n+ i — max{» r : (r,v r ) G 5, r > f n } . (32) 

We compose the envelope out of points densely sampled from 
the line segment connecting (r' n , v' n ) and (r' n+1 ,v' n+1 ), where 
the prime indicates that the coordinates are slightly shifted 
outwards, i.e. r' n — (1 + S)r n and v' n — (1 + S)v n with typical 
S — 0.01. The loop eventually stops after N iterations at 
tn — rup and is followed by an analogous procedure running 
from f o to r\ to complete the envelope along the left part of 
the SOS. The points on the envelope are used as additional 
seeds for the Voronoi tessellation. 

As Fig.|5|shows, the apposition of the boundary points 
as described above ensures that all orbital Voronoi cells are 
closed and confined to an area roughly bounded by the out- 
ermost invariant curve of the SOS. The definition of the 
boundary is purely geometrical and insensitive to numerical 
uncertainties in the orbit integration. The spiky cells along 
the upper boundary belong to seeds of the envelope and do 
not affect the calculation of the orbital phase- volumes. 



3 Note that the Poincare map of the SOS on to itself is area- 
preserving and one would like to have AA S independent of s. The 
Voronoi tessellation however yields only approximately constant 
AA". Nevertheless, as Sec. 151 shows, the resulting phase- volumes 
are of high accuracy. 
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Figure 2. A Voronoi tessellation of the SOS of Fig. Q Open 
circles mark individual intersections of orbits with the SOS, solid 
dots are points added to get the Voronoi cells well behaved at the 
boundaries. 



The Voronoi tessellation used to approximate the inte- 
gral expression in equation 1261 via equation 1291 defines 
a robust method to calculate the relative phase-volume of 
any orbit inside a particular sequence of orbits with com- 
mon E and L z , including resonances and chaotic orbits. 
The areas assigned to the individual orbital imprints in the 
SOS completely fill the area below the ZVC of the SOS. 
Thus, a cruder sampling of the SOS is compensated by 
larger individual orbital phase- volumes. In the limit of an in- 
finitely dense sampling, the assigned "phase-space-weights" 
obtained by the tessellation approach single-orbit phase- 
volumes. 

Calculating AE AL Z . For a complete determination of 
the phase- volumes we also need the relative contribution of a 
whole sequence of orbits with common (E, L z ) as compared 
to other sequences with different energies and angular mo- 
menta. These are described by the factors AL Z AE of the 
orbital phase- volumes (cf. equation 1261 ). They are in fact 
equal for all orbits in the same sequence and need to be 
calculated only once for each sequence. 

Fig. |3| shows an example of the (E, L z )-plane of a typ- 
ical library. The dots display the grid of sampled orbital 
energies and angular momenta. Each dot represents a se- 
quence of orbits with common E and L z but different 73. To 
calculate AL Z AE for a particular sequence, we construct a 
quadrangle around the sequence's (E,L Z ) and estimate the 
product AL Z AE as the surface area enclosed by this quad- 
rangle. The thin lines in Fig.[!|]sliow the boundaries of these 
quadrangles, which are constructed as follows. 

As described in Sec. 12. li the grid in (E, i z )-space is de- 
rived from the requirement that for every pair n < 7j of 
equatorial grid bins, the library contains at least one equa- 



torial orbit with r„ 



n and r a 



a sequence of orbits with (E Ecq ,L z 



- rj. Consider now 
and corresponding 
7" P cri,seq and r ap0iSC q of the equatorial orbit. In (_E,L z )-space 
all sequences inside the boundary of the sampled area are 
surrounded by four other sequences having both their peri- 
and their apocenter in adjacent spatial bins. Let r per ij and 
>"apo, j , 1 ^ J ^ 4, denote the corresponding peri- and apoc- 
enters of the equatorial orbits of these sequences. We con- 



Mapping stationary axisymmetric phase-space distribution functions by orbit libraries 7 







10 



11 



In |E| [km /s 



Figure 3. Typical distribution of sampled energies E and angu- 
lar momenta L z of an orbit library. The thin red lines show the 
boundaries of small cells assigned to each sequence. Their surface 
area is taken to estimate AE AL Z . The potential equals that of 
Figs. Q] and El 



struct a quadrangle around the sequence (E scq ,L z 



by 



connecting the energies and angular momenta of four ficti- 
cious orbit sequences characterised by the pericenter r per i,j 
and apocenter r ap0: j of their equatorial orbits 



^"pcri 



'Deri, i T" ^*d 



peri, j — o V P cri iJ ' pcri,scq 



(33) 
(34) 



and 

- 1, . 

The sequences with the largest apocenters and the 
smallest pericenters, respectively, are surrounded by less 
than four sequences having both their peri- and their apoc- 
enter in adjacent spatial bins. For these sequences we cal- 
culate the edges of the quadrangle as if there were further 
sequences around, whose energies and angular momenta fol- 
low from our spatial grid at smaller radii than r m i n and 
larger radii than r max . 



Sequences with r pC ri,; 



(lying on the upper 



boundary of the sampled area in Fig. [3] and usually con- 
taining only one, approximately circular, orbit) are also not 
surrounded by four sequences as described above. For these 
sequences we take (E seq ,L Z:Scq ) of the actual sequence as 
the upper right edge of the quadrangle. 

As can be seen in Fig. [3] the quadrangles around the 
sequences' energies and angular momenta completely cover 
the sampled part of the allowed area in (E, L z )-space below 
the curve L Z (E) = Lz,cnc- They give a reasonable measure of 
the fractional area in (E, L z )-space, occupied by each orbit 
sequence. 



5 MAPPING DISTRIBUTION FUNCTIONS 
ON TO THE LIBRARY 

In this Section we describe how to use the phase-volumes 
from the previous Section to calculate internal and projected 
properties of stationary DFs using an orbit library. To this 
end, starting with a density profile p and a stationary distri- 
bution function f p connected to p via p — J f p d 3 v, a library 



is constructed as described in Sec. [5] Instead of fitting the 
library to the kinematics of f p , 



LOSVD f (v, os 



fpd v± 



(35) 



we assign an appropriate weight to each orbit such that the 
superposition of all orbits represents f p (see Sec. 13. 21 above'). 
We then compare the internal density distribution put, and 
the anisotropy profile /3nb, as well as the projected kinemat- 
ics GHPnb obtained from the library with the same prop- 
erties p, (3 and GHP directly calculated from the DF (see 
Sec. 15. 11 - 15. 21 . Thereby we can check to which accuracy the 
orbit library reproduces a given dynamical system. 

5.1 Spherical 7- models 

As a first reference case, we explore spherical 7-models. 

Properties of the input model. The stell ar body of the 
reference model is constructed from 7-models ( Dehnen 1993) 
with density 

M r s (3-7) 



(h(r) = 



(36) 



47r ri(r s + r) 4 ~i ' 

They approximate the de Vaucouleurs law of ellipticals 
quite well for 7 6 [1,2]. The DF is assumed to be of the 
Osipkov-Merritt type /q M = foM JE - L 2 /2r 2 a ) JOsipkovl 
ll979UMerrittll985al:lMerrittll985bTl . The corresponding sys- 
tems are isotropic at r < r a and radially anisotropic at 
r » r a : 



0i 



2 1 2 

Vj> + V<t, 

2v 2 



(37) 



We tested various combinations of the parameters (7, r s , r a ). 
However, since the conclusions drawn from the comparisons 
do not depend strongly on 7, the following contains only 
a discussion of the results for the Hernquist model (7 = 
1), where the DF c an be written in terms of elementary 
functions and reads (IHernauistlll99Cl) 



f(E,L)cc 



7 2\5/2 



3arcsing + (1 — 2q 



ly/l^- 



l 4 -8q 2 



3) 4" 



(38) 



with q = ^/r s (E-L 2 /2r 2 a )/GM. 



Comparison of model and library. Fig. 2] shows the 
GHPs, density and anisotropy profiles of a library with ~ 
2 x 8800 orbits, extending from S3 5 x 10" 4 r s to w 28 r s . For 
this library we have used a closed meshed sampling contain- 
ing two different pericenters for each radial bin. The weights 
for the orbits were directly derived from equation 123H and 
the Hernquist DF of equation 138H with r s = 10.5 kpc, a 
total mass of M = 7.5 x 10 11 Mq and r a — 00 (isotropic 
model). The big dots show the expected kinematics, den- 
sity and anisotropy of the Hernquist model. The GHPs were 
obtained by first calculating the LOSVDs at the position of 
the correspo nding spatial bin from the DF using the me thod 
described in lCarollo. de Zeeuw fc van der Mare] 1119951) and 
then fitting a GH-series to the LOSVDs. For the density dis- 
tribution and anisotropy we used equations I36H and 13711 . 
respectively. 
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As the figure shows, the library is able to reproduce the 
GHP and internal density distribution of the model to a high 
degree of accuracy. The mean fractional difference in a is be- 
low Act < 1 per cent and the mean difference in Ha is below 
AH4 < 0.01. The largest deviations between model and li- 
brary occur in the anisotropy profile with rms(/3) = 0.06 
(taken over a whole angular ray). The individual differences 
however are smaller than A/3 = 0.1 over almost the whole 
spatial range covered by the library. Near the inner and outer 
boundary of the library the orbit sampling becomes incom- 
plete with mostly radial orbits coming either from outside 
the outer boundary or from inside the inner boundary are 
missing. Consequently, the library shows a decrease in radial 
velocity dispersion (/3 < 0) as compared to the expectations 
of the isotropic reference model. 

Fig. 13 shows the same as Fig. [I] but for an anisotropic 
Hernquist model with r a = 4 r s . It confirms the results from 
the isotropic model. The offset in the //4-profiles at large 
radii is due to errors in the GH-fit. At these radii the reso- 
lution of our LOSVD-bins is too low to give reliable GHPs. 
However the match of the individual LOSVDs itself is as 
good as at smaller radii. 

Again the largest deviations show up in the /3-profiles, 
with a mean rms(/?hern — Ait>) = 0.03. As in the isotropic 
case the differences between model and library increase when 
approaching the edges of the library, where the radial ve- 
locity dispersion of the library is systematically lower than 
expected. 

5.2 Flattened Plummer model 

We now go one step further and use a flattened test object , 
namely the flattened Plummer model of lLvnden-B ell (1962) 
(normalized such that in the spherical limit M defines the 
total mass of the model) 

AaM) = — ^--[(30 



^)( 



Ua? 



r + a 

22 2, 



+ 



b ) b r cos ($)], 



(39) 



A — (r 2 + a 2 J — 2b 2 r 2 cos 2 (i9). The parameters a and b 
describe the extension of the core and the flattening. The 
part of the distribution function, which is even in L z is 



f p i(E,L z 



V2 
4tt 3 / 2 



T ® DEi + WgLcLlE* 

" ' 15 \ z 



P(- 



(40) 



C oc (3a 2 -26 2 ) and D oc 5b 2 (2a 2 -b 2 ). The Plummer models 
do not rotate as long as the uneven part of / p i vanishes and 
prograde and retrograde orbits exactly balance each other. 

Comparison of model and library. Fig. |S] shows a flat- 
tened Plummer model with a = 5.0 kpc, b — a/2 and M = 
7.5 x 10 1 Mq (big dots) and profiles obtained from a library 
with « 2 x 4400 orbits, extending from « 10 -3 a to « 20 a 
(solid and dashed lines as in Figs.^JandQjJ- The weights were 
derived from / p i via equation 1231 . The kinematics along 
the major and the minor axis have been calculated from 
higher order Jeans equations iMagorrian fc Binnevl ll994'). 
Before determining the GH-parameters the projected mo- 
ments were integrated along a 3.6 arcsec wide major-axis 
slit and a 2.0 arcsec wide minor-axis slit. (Note that for the 
axisymmetric case we take (3 = 1 — a 2 , /a 2 .) 



^ 4 

» 2 

6 
-2 
> -4 



< 



,— , 200 ~- 



6 150 
b 100 

0.04 

0.02 



0.02 

-0.04 

0.1 





-0.5 

log(r/r ) 



2 


r- 1 


III | 1 1 1 1 1 1 1 




fl „ 
q. -2 

M _4 

^ 
-6 




^^•^ 


2 

'"a 
E „ 

60 -4 

* 

-6 




,M_,M, :J:: ,^ 




0.2 



-0.2 


• 










0.2 



-0.2 


• 


• _•-•- -•_•_ • •-*_•_ _• .•_•_ 


• * -•-•-• • • — 



-2 -1 

l°g(r/r„) 



Figure 4. Comparison of a Hernquist model (big dots) and a 
library with weights directly derived from the spherical, isotropic 
Hernquist DF (lines). The upper panel shows the projected kine- 
matics along the major axis (solid line) and minor axis (dashed 
line). The lower panel shows the density distribution (upper two 
rows, [p] = Mq/pc 3 ) and the anisotropy parameter (lower two 
rows) for the minor and major axis, respectively. 



As in the spherical case the Gauss-Hermite parameters 
of the projected kinematics are reproduced to better than a 
few percent. Deviations in the outer parts of the ^-profile 
stem from the GH-fit and are not seen in the LOSVDs. The 
density distribution is also well reproduced up to ~ a/10 
and the anisotropy parameter is \/3\ < 0.1 from the outer 
edge of the library down to as a/10. 
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Figure 5. Same as Fig.|3]but for an anisotropic Hernquist model 
with r a = 4r». 



5.3 Changing the spatial coverage of the library 

The library only discretely represents a finite part of the 
available phase-space. To check how this affects the accuracy 
of the calculation of phasespace integrals of a given DF with 
the library, we did the profile comparisons described in Sects. 
I5.1l and l5.2l for libraries with different spatial extent and for 
different resolutions in the space of orbital integrals. 

The upper panel of Fig. Q shows a and H4 along the 
major and the minor axis for the isotropic Hernquist model 
(big dots). The four different lines show the outcome of four 
libraries with different spatial coverage. For the solid line 
(V m m = 2.5 x 10~ 4 ,r max = 10) (in units of the effective 
radius), for the dotted line (r m i n = 2.5 x 10~ 3 ,r max = 10), 



Figure 6. Comparison of a flattened Plummer model (big dots) 
and a library with weights directly derived from the DF (lines). 
Upper panel: projected kinematics along major axis (solid lines) 
and minor axis (dashed lines). Only moments independent from 
the uneven part of the DF are shown. Lower panel: density dis- 
tribution (upper two rows, [p] = Mqpc -3 ) and the anisotropy 
parameter (lower two rows) for the two axes. 



2.5 X 10" 



5) and 



for the short dashed line (r m i n = 2.0 x lu ",r max 
for the long dashed line (r m i n = 2.5 x 10 _3 ,r max = 5). 

As expected the less extended libraries fail to reproduce 
the innnermost or outermost datapoins, respectively. In the 
vicinity of the equatorial plane (e.g. along the major axis and 
at the central parts of the minor axis) the library becomes 
dominated by azimuthal motion, when approaching r m i n or 
r max , since orbits coming from further outside or inside are 
missing. Consequently, the LOSVDs are too flat (too small 
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Hi) as compared to the expectations (see e.g. the outermost 
parts of the dashed lines of the libraries with small r max 
along the major axis and the innermost parts of the long- 
dashed and pointed lines of the libraries with large r m i n in 
the minor-axis /^-profile) . 

The effect can also be seen in the internal dynamical 
structure, which is illustrated in the lower panel of Fig. Q 
where the anisotropy of the library with respect to <p and i? 
is plotted separately, 



P<f 



1 



<Jr 



1 _ Z± 

or 



(41) 



Near the centre j3 v < along the major and minor axis, 
confirming the dominance of y>-motion brought about by the 
dominance of orbits having their inner turning points there 
and consequently rotate fastly around the axis of symmetry. 
The effect is less pronounced at the outer points of the major 
axis, where the effective potential of the meridional plane 
motion is less dominated by the L 2 -term. 

The /3,9-profiles lack from boundary effects because they 
are independent from the (i?,L z )-sampling and simply re- 
flect the degree to which the SOSs are filled with orbital 
invariant curves. 

Along the minor axis the agreement of library and 
model in projected a is quite good. Near the centre the 
library's a is enhanced due to the orbits having their peri- 
center there. 



5.4 Changing the number of orbits in the library 

Fig.|S]shows the same comparison as Fig.Q but for libraries, 
where we have skipped every second r pcr i resulting in only 
~ 2 x 4700 orbits per library. The gross appearance of Fig. |H] 
is quite similar to Fig.Qwith some minor differences. First, 
the scatter in the GHPs has increased a bit, however the 
match of predictions and library is still on a level of a few 
percent. 

The most striking difference is the increase of radial 
relative to azimuthal motion near the centre of the library. 
Most probably this reflects the fact that the pericenters of 
the orbit sequences are located at the inner edge of each 
radial bin. Therefore the most radial orbits which con- 
tribute also significantly to v^ near their turning points 
move through the whole bin before turning around and thus 
rise the radial velocity dispersion. This effect is strongest in 
the centre since our binning there becomes relatively large 
compared to the variation of the potential. The balance be- 
tween radial and meridional motion is not affected by this 
resolution-effect, because the sampling inside each sequence 
(in the SOS) is independent from the (E, L z )-grid and thus 
independent from the resolution of the sampled peri- and 
apocenters. 



6 FITTING THE LIBRARY 

So far we have omitted the problem of finding the or- 
bital weights Wi according to some given kinematical con- 
straints. This section contains a brief description of our use 
of the maximum entropy technique of Richst one fc Tremaind 
(1988) to fit the library to some LOSVDs. 
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Figure 7. <r and Ha (upper panel) and anisotropy (lower panel) 
along major and minor axis for the isotropic Hernquist model 
(big dots) and four libraries with different spatial extension 
(rminAcff , r max /r eff ): (2.5 X 10~ 4 , 10) solid line, (2.5 X 10~ 3 , 10) 
dotted line, (2.5 X 10~ 4 , 5) short dashed line, (2.5 X 10 -3 , 5) long 
dashed line. 



6.1 Maximum entropy technique 

Given a set of kinematic constraints, we seek the orbital 
weights, that best fit the library to the constraints. These 
weights a re derived from the maximizat ion of an entropy-like 
quantity (Richs tone fc TremaineHl988l) 



S = S-a X 
where 



(42) 
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Figure 8. Same profiles as in Fig.|3 but the libraries have been 
set up with a coarser sampling with roughly half the number of 
orbits as compared to Fig, 171 



2 y^ / LOSVD(lib) 
i,k ^ 



jk 



LOSVD(data) 



jk 



ALOSVD(data) 



jk 



(43) 



gives the departure between the predicted kinematics of the 
library LOS VD (lib) (cf. equation JUJ) and the data kinemat- 
ics LOS VD (data) . Note that the luminosity density v is not 
fitted, but used as a boundary condition (see Richstone et 
al. (in preparation) for details). S is an approximation to 
the usual Boltzmann-entropy 



S= I / lib ln(/ lib )d 3 rd 3 « = 



Wi In 



m 



(44) 



In the absence of any other condition the maximization 
of S enforces the weights Wi to be proportional to the phase- 
volumes Vi. This fact can be used to bias the library towards 
any set of predefined weights. If, for example, we substitute 
the phase- volumes in equation 1441 by Vi — > fi Vi, 



S^S' 



2_, Wi m 



fiVi 



(45) 



then the maximization of S' yields weights Wi proportional 
to fi Vi. According to equation i12,'-j> one can choose the fac- 
tors fi to bias the library towards any given DF /. The 
Boltzmann entropy corresponds to the case of equal a priori 
probabilities fi — fj in phase-space. 

6.2 The smoothing parameter a 

The smoothing parameter a controls the influence of the 
entropy S on the fitted weights. If a is small, then the max- 
imum of 5* is less affected by \ an d the library gives a poor 
fit to the data. Consequently, it will not represent the true 
structure of the object to which it is fitted. If on the other 
hand a is large, then the maximum of S is mostly deter- 
mined by the minimum of y? . In this case the library fits 
the noise in the data. The DF of the library is then highly 
unsmooth and again does not represent the true DF of the 
corresponding object. 

The question of how much smoothing has to be applied 
in order to get an optimal estimate of the true underlying DF 
for a given set of observational data with specific errors and 
spatial sampling will be the content of a forthcoming paper. 
Here, we focus on illustrating the accuracy of our method 
to setup the orbit libraries. In the following we will always 
show the results for that a which gives the best match to 
the input DF. 



7 RECONSTRUCTING DISTRIBUTION 
FUNCTIONS FROM FITTED LIBRARIES 

In this Section we use the DFs of Sec. |3J but instead of 
exploiting equation (12311 to assign the orbital weights and 
to compare spatial profiles of the library and the original 
DF, we now fit the library to the DF as follows. First, we 
calculate the density profile and GHPs connected with the 
DF 



P = 



and 



fd 3 v 



(46) 



LOSVD y (i; los ) = i //dV (47) 

where the GHPs are obtained from the LOSVDs as described 
in Sec. 15. II We compose a library as dscribed in Sec.|2]and fit 
it to the GHPs via the maximum entropy technique of Sec.[£] 
Finally, we compare the orbital weights Wi (a) resulting from 
the fit with those expected from the DF via equation \Y2'M . 
By showing that the fitted weights approximate the input 
DF over a large region in phase-space, we justify that we can 
use the degree to which the library approximates the DF as 
a criterium to determine the optimal amount of smoothing, 
which we will exploit in a subsequent paper in more detail. 
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In order to find the best fit weights that minimize the 
X 2 of Eq, l43L we derived error bars for the LOSVDs by 
first assigning error bars to the GHPs and then determining 
LOSVD errors by means of Monte-Carlo simulations. The 
error for a was choosen to linearily increase with r from 2 
per cent at the innermost data point to 10 per cent for the 
outermost data point. For Hz and H4 the errors increase 
from 0.01 to 0.05. The definition of the errors is somewhat 
arbitrary since we do not add noise to the data points, but 
they are roughly comparable to real data error bars. Since 
v — in the models, the error for v is set to Av — 2 kms^ 1 . 
A detailed investigation of the influence of realistic errors 
on the accuracy of the reconstructed internal properties of 
a fitted library will be presented in a forthcoming paper. 



7.1 Hernquist model 

Fig. |5| shows a comparison of characteristic properties of a 
library fitted to the kinematics corresponding to the dots in 
the upper panel and the original DF. The definition of the 
lines and dots as well as the input DF are the same as for 
Fig. 0] the fit was obtained with a = 0.0046. As expected 
the match to the kinematics and the internal density profile 
is excellent after the fit. The anisotropy is smaller than \/3\ < 
0.1 over a spatial region largely exceeding the area where the 
LOSVDs were fitted. Only near the very centre, the minor- 
axis /3-profile drops significantly because of the lack of radial 
orbits coming from inside the inner boundary of the library 
(cf. Sec. U>J . 

Fig. 1101 shows the DF reconstructed from the fitted 
weights via equation (1221 (dots) together with the input 
DF (thick line). Each dot represents the phase-space den- 
sity along one single orbit, the densities are scaled according 
to y^ Wj — y^ Vj — 1 . Over a region covering 90 per cent 
of the library's mass, the rms-difference between the Hern- 
quist DF and the orbital phase-space densities is 12.1 per 
cent. The remaining departures between model and fit are 
mostly due to boundary effects arising from the discrete and 
finite nature of the library. 

Fig. 1111 shows the fractional differences of input model 
and library as a function of orbital energy E and angular 
momentum L z . For each dot, the contributions of individual 
orbits with common E and L z have been integrated. Larger 
dots correspond to larger differences between input DF and 
fitted library. From Fig. Illl one sees that the remaining de- 
viations between library and input DF mostly stem from 
orbits lying at the boundary of the phase-space-region cov- 
ered by the library. Since the library only contains a finite 
number of all orbits, the fit to the kinematics with the den- 
sity as a boundary condition enforces some redistribution 
of orbits as compared to the original DF. For example, at 
the outer boundary of the library (E « 0) the fitted orbital 
phase-space densities are too large as compared to the in- 
put DF. These orbits compensate the cut-off in energy and 
contain all the light that should have been distributed along 
even lower bound orbits. For the same reason, the library 
fails to reproduce the Hernquist DF near the most bound 
orbits. 

Fig. H2l shows the results when fitting the same library 
to the projected kinematics of the anisotropic Hernquist 
model with r a — 4r s , corresponding to the dots in the up- 
per panel of the Figure. Again, after the fit the library per- 
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Figure 9. Comparison of a library fitted to the LOSVDs of a 
spherical, isotropic Hernquist DF (lines) and the Hernquist model 
itself (big dots). The upper panel shows the projected kinematics 
along the major axis (solid line) and minor axis (dashed line). 
The lower panel shows the density distribution (upper two rows, 
[p] = Mq pc -3 ) and the anisotropy parameter (lower two rows) 
for the two axes. 



fectly reproduces the internal density profile and the pro- 
jected kinematics. The mismatch in the outer parts of the 
//4-profiles result from errors in the GHP-fit (cf. Sec. I5.1|l . 
However, we don't fit the library to the GHP, but directly to 
the LOSVD. The /3-profiles of the library follow the expected 
curves well inside the region covered with kinematical con- 
straints. In the outer parts however they do not follow the 
input model to predominant radial motion but turn back to 
an isotropic appearance. This is a reflection of the entropy 
maximization used in the fit, which forces those parts of 
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-E [kmVs 2 ] 

Figure 10. Comparison of the DF of a spherical Hernquist model 
(solid line, units defined in the text) with the phase-space densi- 
ties obtained from a library fitted to GHPs along two perpendic- 
ular axes in the galaxy (details in the text) . Each dot represents a 
single orbit. The rms between library and model is 12.1 per cent 
over a region covering 90 per cent of the library's mass. 
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Figure 11. The fractional difference between a spherical Hern- 
quist model and a fitted library as a function of orbital energy E 
and z-angular momentum L z . Larger dots corresponds to larger 
differences. For open dots, the DF of the library overestimates 
the real DF, for the solid dots it underestimates the DF. 



the library to isotropy, which are not constrainted by data 
points. 

To confirm that, we refitted the library, however re- 
placed the Vi in equation 1)440 by the weights of the 
anisotropic Hernquist DF following from equation l)23|l . 
Since for the maximum entropy solution of equation 1441 
(without any other condition) the weights lOj are pro- 
portional to the values Vi, now being the weights of the 
anisotropic DF instead of the phase- volumes, this biases the 
fit towards the anisotropic Hernquist model. The character- 
istics of the corresponding fit are displayed by the dotted 
lines in Fig. 1121 The projected kinematics and internal den- 
sity are indistinguishable from the maximum entropy fit, but 
now the anisotropy profile is in perfect agreement with the 
input model. 




Figure 12. Same as Fig.|[5]but for an anisotropic Hernquist model 
with r a = Ar s . The dotted line shows the result of a fit with 
"biased weights" (see text for details). 



7.2 Flattened Plummer model 

Fig. llSl shows the GHPs and internal density and anisotropy 
of the Plummer model with b = a/2 of Sec. l5.2l together with 
a fitted library containing ~ 2 x 4400 orbits. The library was 
fitted to the LOSVDs corresponding to the dots of the upper 
panel of the Figure with a smoothing parameter of a w 0.03. 
The small deviations between the library's kinematics and 
the model in the upper panel of the figure are due to the low 
resolution in the GH-fit and are not seen in the LOSVDs 
used for the fit. The anisotropy parameter is confined to 
\P\ < 0.1 all over the region where the library is constrained 
by kinematic data. 
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Figure 14. Deviations from the reconstructed DF of a fitted 
library and the Plummer DF of Fig. 1131 Each dot represents one 
sequence of orbits with common E and L z . For the open dots, the 
DF of the library overestimates the real DF, for the solid dots it 
underestimates it. Larger dots indicate larger differences. 




-6 

0.4 

0.2 



-0.2 

-0.4 

0.4 

0.2 



-0.2 

-0.4 







• — • • » • 



-H- 



-1 0.5 0.5 1 

log(r/a) 



Figure 13. Comparison of a flattened Plummer model (big dots) 
and a fitted library (lines). The upper panel shows the projected 
kinematics along the major axis (solid lines) and minor axis 
(dashed lines). The lower panel shows internal moments along 
the minor and major axis, respectively (units as in Fig. Itjl . 



The rms-difference between the reconstructed DF and 
the input model is fa 15 per cent over a region covering 90 
per cent of the library's mass. As Fig. H4l shows. differences 
between model and library are confined to the boundaries 
of the sampled (E, L z )-region of the phase-space. As for the 
Hernquist model, the reason for these differences is the in- 
complete orbit sampling at the edges of the library. 



8 SUMMARY 

We have presented a modified version of the Schwarzschild 
code of Richstone et al. (in preparation). The code involves 
a new orbit sampling at given energy E and angular mo- 
mentum L z and a new implementation for the calculation 
of the orbital phase- volumes. 

For our libraries we supplement the drop of orbits with 
common energy E and angular momentum L z from the ZVC 
as described in Richstone et al. (in preparation) by scanning 
the SOS with a resolution that varies as the sampling pro- 
gresses from the more radial to the more shell-type orbits. 
This sampling has been shown to completely fill the SOS 
connected with a pair (E,L Z ) with orbital imprints. 

A Voronoi tessellation of the SOSs of orbits with com- 
mon E and L z allows us to calculate the phase-space- 
volumes of individual orbits in any axisymmetric potential. 
With the phase- volumes we can convert the orbital weights 
describing the relative contribution of the orbits to the whole 
library into phase-space densities and vice versa. As a first 
application we use the densities to check our method to setup 
the library in two different ways. 

First, we calculate spatial profiles of internal and pro- 
jected properties of isotropic and anisotropic DFs of spher- 
ical 7-models as well as of the flattened Plummer model 
with the library. The density profiles, anisotropy profiles and 
projected kinematics of the library closely match those di- 
rectly inferred from the corresponding DF. The errors in 
the higher order GH-parameters H n are AH n < 0.01, for 
n = 3, 4 and the fractional error in the projected dispersion 
Act < 1 per cent are accurate on a level better than that of 
present day observational errors. The largest deviations oc- 
cur in the anisotropy profile, but are smaller than A/3 < 0.1 
at almost all positions in the library, however increase to- 
wards the edges of the spatial region which is covered by 
the orbits. This boundary effect is caused by the locally in- 
complete orbit sampling there. If in practical applications 
the libraries are constructed to extend beyond the area with 
observational constraints, then these inaccuracies are negli- 
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gible. Hence, our libraries fairly represent the phase-space 
structure of the considered models. 

As a second application we fitted libraries to the GHPs 
of the same spherical 7- models and flattened Plummer mod- 
els. The reconstructed DFs match the input DF with a rms 
of about 15 per cent over a region covering 90 per cent of 
the library's mass. The remaining deviations are mostly re- 
stricted to orbits at the boundary of the phase-space volume 
represented by the library. This is not unexpected since the 
library only discretely represents a finite subregion of the 
input system. Consequently some redistribution of orbits is 
necessary to compensate for orbits not included in the li- 
brary. 

We will investigate the influence of observational errors 
on the reconstructed DFs and of the amount of smoothing 
applied in the fit in a forthcoming publication. In a further 
step we will reconstruct the internal structure and mass com- 
position of a sample of flattened early-type galaxies in the 
axisymmetric approximation. 
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